RaDAr

In this section we will use the library RadAR (Radiomics Analysis with R). This is a software to perform comprehensive analysis of radiomic features. RadAR allows the users to carry out the entire processing of radiomic datasets, from data import to feature processing and visualization and implements multiple statistical methods for the analysis of these data.

First we load the RadAR library.

library(RadAR)

With RadAR we can import the radiomic features that we extracted before using pyradiomics.

TO DO: explain normalizacion y shift, y con normalizacion.

We will create the rdr object with normalization and shift, and another one with normalization only. We use the function import_pyradiomics(path_to_folder) to import the results we obtained using pyradiomics, where path_to_folder is the path to the folder containing the pyradiomics csv files. The function filters out duplicated feature names.

rdr_NormShift <- import_pyradiomics(dir="/Users/andrealetaalfonso/Desktop/TFG/A_Project/results/ECLIPSE_results_ParamsNormShiftCSV")
rdr_Norm <- import_pyradiomics(dir="/Users/andrealetaalfonso/Desktop/TFG/A_Project/results/ECLIPSE_results_ParamsNormCSV")
rdr_NormShift
## class: SummarizedExperiment 
## dim: 1232 2580 
## metadata(1): extractor
## assays(1): values
## rownames(1232): Elongation.original Flatness.original ...
##   Complexity.wavelet.LLL Strength.wavelet.LLL
## rowData names(4): feature_name image_type feature_description
##   feature_type
## colnames(2580): 023140000001_INSP_STD_L1_ECLIPSE.nii.gz.tsv
##   023140000002_INSP_STD_L1_ECLIPSE.nii.gz.tsv ...
##   026658006357_INSP_STD_L3_ECLIPSE.nii.gz.tsv
##   026658006358_INSP_STD_L2_ECLIPSE.nii.gz.tsv
## colData names(3): filename sample_id mask_id
rdr_Norm
## class: SummarizedExperiment 
## dim: 1232 2580 
## metadata(1): extractor
## assays(1): values
## rownames(1232): Elongation.original Flatness.original ...
##   Complexity.wavelet.LLL Strength.wavelet.LLL
## rowData names(4): feature_name image_type feature_description
##   feature_type
## colnames(2580): 023140000001_INSP_STD_L1_ECLIPSE.nii.gz.tsv
##   023140000002_INSP_STD_L1_ECLIPSE.nii.gz.tsv ...
##   026658006357_INSP_STD_L3_ECLIPSE.nii.gz.tsv
##   026658006358_INSP_STD_L2_ECLIPSE.nii.gz.tsv
## colData names(3): filename sample_id mask_id

The rdr object is a SummarizedExperiment.

TO DO: explain more & add image SummarizedExperiment

The available clinical data can be also imported in the rdr object.

We can extract the data with the function assay().

data <- assay(rdr_NormShift)
data[1:10,1:2]
##                                  023140000001_INSP_STD_L1_ECLIPSE.nii.gz.tsv
## Elongation.original                                             7.477679e-01
## Flatness.original                                               6.687235e-01
## LeastAxisLength.original                                        1.879875e+02
## MajorAxisLength.original                                        2.811139e+02
## Maximum2DDiameterColumn.original                                2.941598e+02
## Maximum2DDiameterRow.original                                   2.479516e+02
## Maximum2DDiameterSlice.original                                 2.735781e+02
## Maximum3DDiameter.original                                      2.943569e+02
## MeshVolume.original                                             5.047730e+06
## MinorAxisLength.original                                        2.102080e+02
##                                  023140000002_INSP_STD_L1_ECLIPSE.nii.gz.tsv
## Elongation.original                                             7.427334e-01
## Flatness.original                                               5.885148e-01
## LeastAxisLength.original                                        1.577146e+02
## MajorAxisLength.original                                        2.679874e+02
## Maximum2DDiameterColumn.original                                2.733002e+02
## Maximum2DDiameterRow.original                                   2.552117e+02
## Maximum2DDiameterSlice.original                                 2.520496e+02
## Maximum3DDiameter.original                                      2.797088e+02
## MeshVolume.original                                             3.359750e+06
## MinorAxisLength.original                                        1.990432e+02
data2 <- assay(rdr_Norm)
data2[1:10,1:2]
##                                  023140000001_INSP_STD_L1_ECLIPSE.nii.gz.tsv
## Elongation.original                                             7.477679e-01
## Flatness.original                                               6.687235e-01
## LeastAxisLength.original                                        1.879875e+02
## MajorAxisLength.original                                        2.811139e+02
## Maximum2DDiameterColumn.original                                2.941598e+02
## Maximum2DDiameterRow.original                                   2.479516e+02
## Maximum2DDiameterSlice.original                                 2.735781e+02
## Maximum3DDiameter.original                                      2.943569e+02
## MeshVolume.original                                             5.047730e+06
## MinorAxisLength.original                                        2.102080e+02
##                                  023140000002_INSP_STD_L1_ECLIPSE.nii.gz.tsv
## Elongation.original                                             7.427334e-01
## Flatness.original                                               5.885148e-01
## LeastAxisLength.original                                        1.577146e+02
## MajorAxisLength.original                                        2.679874e+02
## Maximum2DDiameterColumn.original                                2.733002e+02
## Maximum2DDiameterRow.original                                   2.552117e+02
## Maximum2DDiameterSlice.original                                 2.520496e+02
## Maximum3DDiameter.original                                      2.797088e+02
## MeshVolume.original                                             3.359750e+06
## MinorAxisLength.original                                        1.990432e+02

We can check which image types (i.e., normal or wavelet filtered) features have been extracted from.

print_image_type(rdr = rdr_NormShift)
## [RadAR] Available image types are
## original
## log.sigma.1.0.mm.3D
## log.sigma.2.0.mm.3D
## log.sigma.3.0.mm.3D
## log.sigma.4.0.mm.3D
## log.sigma.5.0.mm.3D
## wavelet.LLH
## wavelet.LHL
## wavelet.LHH
## wavelet.HLL
## wavelet.HLH
## wavelet.HHL
## wavelet.HHH
## wavelet.LLL

The rdr object can be also filter by feature type.

print_feature_type(rdr = rdr_NormShift)
## [RadAR] Available feature types are
## first_order_shape
## first_order_statistics
## grey_level_co-occurrence_matrix
## grey_level_dependence_matrix
## grey_level_run_length_matrix
## grey_level_size_zone_matrix
## neighbouring_grey_tone_difference_matrix

To access to the elements of patients (columns) or features (rows) we can do:

names(colData(rdr_NormShift)) # columns
## [1] "filename"  "sample_id" "mask_id"
names(rowData(rdr_NormShift)) # rows
## [1] "feature_name"        "image_type"          "feature_description"
## [4] "feature_type"

Correlation matrix

Visualize the correlation matrix:

plot_correlation_matrix(rdr = rdr_NormShift, 
                        method_correlation = "spearman", 
                        view_as = "heatmap", 
                        which_data = "normal")

A different way to visualize the correlation matrix is by using correlation plot:

plot_correlation_matrix(rdr = rdr_NormShift, 
                        method_correlation = "pearson", 
                        view_as = "corrplot", 
                        which_data = "normal")

Pre-processing and visualization of radiomic features

Scale features:

rdr_scaled <- scale_feature_values(rdr = rdr_NormShift, method = "minmax")

Group features:

rdr <- do_hierarchical_clustering(rdr = rdr_scaled, 
                                  which_data = "scaled", 
                                  method_dist_col = "correlation.spearman")

Heat map:

plot_heatmap_hcl(rdr = rdr)

We can also save the rdr object to .rda.

save(rdr_NormShift, file = "rdr_Andrea_NormShift.rda")
save(rdr_Norm, file = "rdr_Andrea_Norm.rda")

We can open the .rda file with the following code:

library(admisc)
## 
## Attaching package: 'admisc'
## The following object is masked from 'package:Matrix':
## 
##     expand
## The following object is masked from 'package:S4Vectors':
## 
##     expand
## The following objects are masked from 'package:igraph':
## 
##     intersection, simplify
test_NormShift <- listRDA("rdr_Andrea_NormShift.rda")
test_NormShift
## $rdr_NormShift
## class: SummarizedExperiment 
## dim: 1232 2580 
## metadata(1): extractor
## assays(1): values
## rownames(1232): Elongation.original Flatness.original ...
##   Complexity.wavelet.LLL Strength.wavelet.LLL
## rowData names(4): feature_name image_type feature_description
##   feature_type
## colnames(2580): 023140000001_INSP_STD_L1_ECLIPSE.nii.gz.tsv
##   023140000002_INSP_STD_L1_ECLIPSE.nii.gz.tsv ...
##   026658006357_INSP_STD_L3_ECLIPSE.nii.gz.tsv
##   026658006358_INSP_STD_L2_ECLIPSE.nii.gz.tsv
## colData names(3): filename sample_id mask_id
test_Norm <- listRDA("rdr_Andrea_Norm.rda")
test_Norm
## $rdr_Norm
## class: SummarizedExperiment 
## dim: 1232 2580 
## metadata(1): extractor
## assays(1): values
## rownames(1232): Elongation.original Flatness.original ...
##   Complexity.wavelet.LLL Strength.wavelet.LLL
## rowData names(4): feature_name image_type feature_description
##   feature_type
## colnames(2580): 023140000001_INSP_STD_L1_ECLIPSE.nii.gz.tsv
##   023140000002_INSP_STD_L1_ECLIPSE.nii.gz.tsv ...
##   026658006357_INSP_STD_L3_ECLIPSE.nii.gz.tsv
##   026658006358_INSP_STD_L2_ECLIPSE.nii.gz.tsv
## colData names(3): filename sample_id mask_id